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Chapter  1 
Introducti  on 

1.1  In  many  physical  problems  of  interest,  with  sonar,  radar,  and  seismology 
as  examples,  the  time  records  of  the  outputs  of  an  array  of  sensors  are  observed 
over  some  time  interval  and  used  to  estimate  the  position  of  a  distant  noise 
source.  Inversely,  the  position  of  the  distant  noise  source  may  be  known,  and  the 
intent  may  be  to  estimate  the  positions  of  the  sensors  comprising  the  array. 
Typically,  the  sensor  outputs  are  amplitude  scaled  and  delayed  replicas  of  the 
waveform  from  the  distant  noise  source,  corrupted  by  additive  noises,  usually 
local  in  origin.  If  the  amplitude  gradient  across  the  array  of  the  waveforms  from 
the  distant  source  is  negligible,  essentially  all  of  the  geometric  information  is 
encoded  in  the  set  of  delays  associated  with  the  propagation  across  the  array  of 
the  wavefronts  from  the  distant  source.  This  paper  discusses  the  theoretical 
bounds  on  the  precision  with  which  the  set  of  delays  can  be  measured,  and  shows 
that  either  properly  filtered  beamformers  or  properly  filtered  systems  of  correla¬ 
tors  can  be  i^ed  to  obtain  estimates  that  achieve  the  theoretical  bound. 

1.2  The  theoretical  bound  discussed  is  the  Cramlr-Rao  matrix  bound  [1], 
which  is  the  appropriate  bound  to  use  when  large  numbers  of  samples,  or 
equivalently,  long  observation  times,  are  used. 

1.3  For  the  purposes  of  this  paper,  it  is  more  convenient  to  use  its  inverse, 
the  Fisher  Information  Matrix,  and  to  compare  the  inverses  of  the  matrices  for  the 
beamformer  and  multiple  correlator  delay  measurement  schemes  to  the  Fisher 
Information  Matrix. 
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1.4  The  results  of  Chapters  2  and  3  are  developed  in  greater  detail  in  [2], 
as  is  the  first  half  of  Chapter  4.  The  remainder  of  Chapter  4  and  Chapter  5  have 

not  been  reported  elsewhere. 

1.5  The  following  notation  is  used:  If  A  is  a  matrix,  A‘^  is  its  inverse, 

ic  T 

A  is  its  conjugate,  A  is  its  transpose,  tr  A  is  its  trace,  and  det  A  is  its 
determinant.  A  square  matrix  whose  elements  off  the  main  diagonal  are  all  zero 
may  be  written  as  diag  (a-j,  a2,  ...  ,  an),  where  a^  is  the  i-th  diagonal  element. 
Vectors  are  column  vectors  unless  otherwise  specified.  ]_  denotes  a  vector  with 
every  element  a  one  (1).  0  is  a  matrix  of  zeros,  and  I  is  the  identity  matrix. 

<•>  is  the  expectation  operator,  and  grad  f  is  the  row  vector  which  is  the  gradient 
of  the  scalar  f.  The  gradient  of  a  vector  is  the  matrix  in  which  the  i-th  row  is 
the  gradient  of  the  i-th  component  of  the  vector.  The  Kronecker  delta  is  denoted 
in  the  usual  way  as  6^-.  Integrals  of  the  form  fly 2  f  dt  are  written  as  Ij  f  dt, 

Q)., 

and  integrals  of  the  form  /_^f  d(u  are  written  as  / g  f  dw.  The  quantity  is 

defined  as  =  Nwg,  where  N  and  uq  are  defined  in  Chapter  2. 

1.6  The  symbols  MLE,  FIM,  CRMB,  and  HOT  are  abbreviations  for  Maximum 
Likelihood  Estimate,  Fisher  Information  Matrix,  Cramer-Rao  Matrix  Bound,  and  Higher 
Order  Terms  (as  in  series  expansions),  respectively.  The  symbols  (CRMB)  and  (FIM) 
represent  the  designated  matrices  themselves. 


1-2 


NOLTR  72-220 


-3- 


t ' 


Chapter  2 

The  Fisher  Information  Matrix 

2.1  Assume  that  the  signal  wave  fronts  from  a  distant  noise  source  propagate 
across  an  M  element  array  of  sensors  and  that  the  signal  amplitude  gradient  across 
the  array  is  negligible.  The  signal  at  the  i-th  sensor  is  s(t-d^),  where  s(t)  is 
the  signal  at  a  reference  point  near  the  array,  and  d^  is  the  delay  at  the  i-th 
sensor.  Without  loss  of  generality,  the  reference  point  is  assumed  to  be  at  the 
location  of  the  first  element  in  the  array.  Thus,  d-j  =  0.  The  output  of  the  i-th 
sensor  is 

xi (t)  =  s(t-di )  y  ni (t) ,  (1) 

where  ( t)  is  the  additive  sensor  noise.  The  M  sensors  are  observed  for  T  seconds, 

-  T/2  <  t  <  T/2,  and  the  M  time  records  are  represented  by  Fourier  coefficients, 

X- (wk}  =~-  fj  x. (t)  exp{-jko>0t]  dt,  (2) 

where  a>0  =  2tt/T,  and  =  ku>0.  The  following  assumptions  will  determine  the  joint 
density  function  for  the  random  Fourier  coefficients. 

a.  The  random  signal  and  each  of  the  M  additive  sensor  noises  are  all 
stationary  zero-mean  Gaussian  random  processes. 

b.  All  of  the  random  processes  are  independent. 

c.  T  is  large  compared  to  the  correlation  times  of  the  random  processes, 
and  also  to  the  time  needed  for  the  signal  wave  fronts  to  transverse  the  array. 

2.2  If  X  is  a  vector  containing  the  Fourier  coefficients  as  elements,  if  S(to) 
and  N-U)  are  the  signal  and  noise  power  spectra  at  the  i-th  sensor,  and  if  only 
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Fourier  coefficients  up  to  frequency  Nu0  are  to  be  processed,  the  density  for  X 
can  be  written  as 


u  n 

p(X)  =  n  det  R(k)]"1  exp{-  l  XT(k)R'1(k)X*(k)}, 


where: 


X(k)  =  {X, (««k).  x2(»k) . X^))1 


V(k)  -  (l.e^  ...  ,ej“kdM)T 
N(k)  =  diag(N^(a)k),  N2(wk),  ...  , 

R(k)  =  K(k)  +  s(»k)  V*(k)  VT(k) 


W> 


2.3  In  what  follows  the  frequency  arguments  of  the  functions  discussed  will 

be  generally  suppressed.  7  will  indicate  a  sum  over  the  positive  Fourier 

Ba¬ 

ft 

frequencies  being  considered,  while  [  will  be  tsed  to  denote  the  array  sum  J  . 

i  i=l 


When  a  appears  following  a  sum  ^  ,  it  will  be  understood  to  stand  for  kwQ. 

2.4  Since  det  R(k)  =  (1  +  £.  S/N - )  det  N(k),  only  the  exponential  part  of 
the  density  function  will  depend  on  the  d^.  Let  the  signal  delay  vector  D  be 


defined  as 


u  -  v«2»  u3» 

2.5  The  likelihood  ^unction  for  0  is 


D-  (d9,  d~,  ...  ,  dM)  . 


L(D)  =  tJ*  n  det  R]"1  c-xp{-  l  xWl. 
B+  B+ 


m 
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2.6  The  CRMB  for  unbiased  estimators  of  the  vector  argument  of  the  likelihood 
function  is  the  inverse  of  the  FIM,  denoted  by  (FIM),  where 

(FIM)  =  -  <  grad(grad  In  L(D))^  >.  (7) 

The  gradients  in  equation  (?)  are  taken  with  respect  to  the  components  of  the 
vector  D.  If  6  is  defined  by 


6  =  S/(l  +  l  S/N.), 


the  inverse  of  the  matrix  R  is 


-l  -l  _i  *  t  «.i  ,  , 

R  1  =  N  1  -  GN  1  V  V1  N  (9) 

Provided  that  only  the  elements  of  0  depend  on  a  and  b,  the  typical  element  of  the 

FIM  has  the  form 

-  <  af  a! *-(D)  >  *  -  l  G  <  xV1  ^  V*  vVY  > 

B+ 

o  c  a{D.-D_)  3(D,-DJ 

i  “Gn"k~5T^“^-  00) 

B+  km  *  m 

2.7  From  equation  (10)  it  follows  that  the  FIM  pertinent  to  the  estimation  of 


the  vector  D  is 


(FIM)  =  [  ?J - 

3+  1  +  I  IT 

1  t 


(tr  N’1)  N*1  -  HpWNp1]. 


In  equation  (11),  N~^is  the  matri  with  the  first  row  and  column  partitioned 
away.  Because  Gf  the  assumed  smoothness  of  all  of  the  spectra  relative  to  the 
frequency  increment  u)Q  =  2ti/T,  the  FIM  can  also  be  written  as 

<«">  *  H  k  “2  7-4t-  [(tr  "■’>  "p1  -  Np’  1 1T  Np’3  «“• 

1  +Iht 
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Chapter  3 

The  Maximum  Likelihood  Estimate 

3.1  It  is  well  known  that  when  the  MLE  is  based  on  a  large  number  of 
independent  samples,  it  is  consistent,  asymptotically  normal,  ard  asymptotically 
efficient  [3].  Since  the  observation  time  T  is  large  compared  to  the  process 
correlation  times,  there  should  be,  in  some  seise,  a  large  number  of  independent 
samples.  The  covariance  matrix  for  the  error  in  the  MLE  for  D  should  be  the  CRMB, 
at  least  to  first  order. 

3.2  The  results  that  follow  are  independent  of  the  true  de^ay  vector,  and 
the  equations  for  the  likelihood  function  and  MLE  are  considerably  simplified  if 
the  true  delay  vector  is  assumed  to  be  0.  The  vector  D  of  this  chapter  is  the  MLE 
and  is  therefore  the  measurement  error.  The  steering  vector  corresponding  to  the 
error  D  is 

VT=  (1,  expijudg),  ...  ,  exp{judM>).  (13) 

The  MLE  vectors  D  and  V  satisfy 

0  =  grad  In  L(D) 

X-X* 

=  grad  \  G  H  ^n-oxp{ju.(dn  -  d^}  (14) 

B+  in  n 

=  grad  (A  +  BD  +  DTCD  +  HOT), 

where  by  expanding  exp{ju(dn  -  d.)}  as  a  power  series,  the  vector  B  is  seen  to  be 

B  =  l  jwGlV^XX*1  -  xVjT1,  (15) 

8+  v  v  v 
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while  the  matrix  C  is  determined  by 

^D^CD  s  ][  ^jw)^6  ££  ITTT^^k  "  •  (16) 

B+  Ik  1  k 

In  equation  (15),  Xp  is  the  single  frequency  data  vector  with  the  first  element 

* 

partitioned  away.  The  terms  X.X^  in  equation  (16)  are  elements  of  the  sample 

covariance  matrix  at  a  single  frequency  based  on  T  seconds  of  data.  These  sample 

covariance  elements  do  not  converge,  even  if  T  is  arbitrarily  long  [4].  However, 

since  T  is  large  compared  to  the  process  correlation  times,  the  spectra  are 

smooth  enough  so  that  the  sample  covariance  can  be  averaged  with  samples  from 

nearby  frequencies  to  provide  statistical  convergence.  The  £  summation  in  equa- 

8+ 

*  * 
tion  (16)  provides  such  an  averaging  of  the  X^X^.  Thus,  it  is  assumed  that  X-X^ 

can  be  replaced  by  R-k  =  <  X..x£  >  in  equation  (16),  from  which  it  follows  that 
C  =  <  C  >,  or 


2rT7T-  CtrN-’lN-’-N-'llV1] 


C  =  -l  2u  - 

t  Wfr 

i  1 


=  -  (FIM). 

3.3  From  equation  (15),  it  follows  immediately  that  <  B  >  =  0,  and  not  so 
immediately  that 


<  BTB*  >  *  l  2u> 


2  S‘ 


B+  Wt 

i  1 


[(tr  N*1)  N"1  -  N"1  llTNp1] 


=  (FIM). 


3.4  Neglecting  the  HOT,  and  assuming  C  =  <  C  >,  the  vector  D  is  given  by 

D  =  -  <  C  >-V  ,  < 
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so  that:  . 

i 

<  !D  >  =  0  ,  t  :  (20) 

and  '  ,  ’  , 

<  dV  >  =  <  C  !>_1  <  BTB*  >*  <  l  >']  t 
'  -  (FIM)*1 

(21) 

=  (CRMB)  .  ' 

3.5  Thus  the  MLE  is  unbiased,  and  in  the  limit  of  large  T  achieves  the  CRMB. 

i  * 

The  MLE  processor  can  be  readily  implemented  as  indicated  in  Figure  1.  The  MLE 

processor  is  just  a  steered  and  filliered  beamformer  followed  by  a  square-law  . 

averager.  The  individual  inputs  are  each  steered  and  then  filtered  with  filters 

whose  frequency  responre  is  the  inverse  of  the  additive  noise  ‘pectrum  for  that 

particular  sensor.  The  beam  sum  is  then  formed  and  fed  to  a  filter  whose  squared 

,  2‘ 

magnitude  response  is  f F( jr^u)  j  =  G(u).  The  MLE  is  determined  as  that  set  of 
steering  delays  that  gives  the  maximum  deflection  of  the  output  meter. 


i  i 


x,  (t) 


STEERING 

DELAYS- 


fIGURE  1:  Beamformer  Implementation  of  the  MLE  Processor 


NOLTR  72-220 


1 


Chapter  4 

Correlator  Delay  Measurement  Systems 

4.1  The  following  scheme  can  be  used  to  estimate  the  unknown  delay  vector  D. 
Let  a  system  of  correlators  be  used  to  form  all  the  M(M-i)/2  correlograms  corre¬ 
sponding  to  processing  all  of  the  M  input  wave  forms  taken  two  at  a  time.  The 
individual  correlators  are  assumed  to  have  Input  filters  for  each  channel*  and  the 
position  of  the  correlogram  peak  is  used  as  a  signal  delay  estimate  for  that 
sensor  pair.  If  each  correlator  is  to  provide  an  unbiased  estimate  of  the  corre¬ 
sponding  signal  delay,  the  input  filters  must  have  the  same  phase  response,  and 
hence  can  be  taken  to  be  identical  filters.  A.  typical  correlator  is  shown  in 
Figure  2.  The  steering  delay  is  adjusted  to  give  the  maximum  deflection  of  the 
meter,  and  this  defines  the  delay  estimate  for  the  correlator. 


STEERING 

DELAY 


FIGURE  2:  A  Typical  Multiplier  Correlator 
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4.2  Let  d..  be  the  correlator  estimate  fo**  the  signal  delay  from  the  i-th 

'  J 

to  the  j-th  sensor,  based  on  the  correlation  of  the  X.(t)  and  X.(t)  time  records. 

Let  e- .  and  F^.  be,  respectively,  the  error  in  the  estimate  d-j,  and  the  filter 
used  on  the  inputs  to  the  correlator.  Define  the  following  scalars,  vectors,  and 
matrices: 


G(ij;kl)  =  C(S+Ni5ik)(S+Nj6jl)  -  (S+N-s^MS+M^)] 
DC  =  ^d12*  d13*  ***  *  dlH*  d23*  ***  *  d(M-l)M^ 


(e12*  e13* 


,  « 


(M-l)M 


)' 


F  =  diag  (jF-^l  *  ^13^  *  ***  *  ^(M-l)M^  ^ 


(22) 


^  ~  Js 


2-c  x 
i;  dr  Cju 


'  diag  (K1J?,  K^,  ...  , 

G  =  [G(ij;kl)] 

4.3  The  square  matrix  G  has  for  its  elements  the  scalars  G(ij;kl)  positioned 
according  to  the  scheme  determined  by  the  order  of  the  subscripts  in  EET,  where  ij 
is  the  row  designation,  and  kl  the  column  designation. 

4.4  Using  equations  (22),  the  covariance  matrix  for  the  correlator  scheme 
measurement  error  vector  can  be  compactly  written  [2]: 


PE  =  <  EET  > 

=  K-%  w^FGFdu)  K'1. 


(23) 


4.5  The  correlator  delay  measurement  vector  is  related  to  the  vector  to  be 
estimated  ,D,  by  the  equation 


Dc  =  AD  +  E, 


(24) 


4-2 


WLTR  72-220 

where  the  matrix  A,  with  its  rows  and  columns  labeled  with  the  same  sets  of  ordered 
subscripts  used  for  the  elements  of  Dj.  (for  the  rows)  and  D  (for  the  columns),  has 
for  its  element  in  the  i  j,k  position 

A(ij;k)  -  6jk  -  6.k.  (25) 

4.6  Since  <  E  >  =  0,  the  Gauss-Markov  estimate  [5]  for  the  vector  D  based  on 
the  correlator  measurements  is 

5  *  CATP"1A3“1ATDc.  (26) 

and  the  covariance  matrix  for  the  Gauss-Markcv  estimate  is 

<  (D-D)(D-D)T  >  =  [A^A]'1.  (27) 

4.7  In  [2]  it  is  demonstrated  that  if  M  =  2  or  M  =  3,  the  Gauss-Markov 
estimate  for  D  achieves  the  CRMB  provided  that  the  filters  satisfy 

iFjjl2  =  s2/NiNj  .  (28) 

1  +  l  ~ 

kN* 

It  was  conjectured  in  [2]  that  the  choice  of  filters  given  by  equation  (28)  was 
optimum  for  M  >  3.  The  conjecture  is  in  fact  valid,  and  the  Gauss-Markov  estimate 
so  obtained  is  in  fact  efficient,  that  is,  achieves  the  CRMB.  This  is  demonstrated 
in  what  follows. 

4.8  From  the  definitions  in  equations  (22),  and  with  the  filters  defined  by 
equation  (28),  the  FIM  of  equation  (11)  is 

(FIM)  =  ATKA  „  (29) 

which  suggests  investigating  the  possibility  that 

aVa  «=  ATKA.  (30) 
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With  the  filters  given  by  equation  (28),  the  matrix  FSF  can  be  written  as 

2 

FSF  *  — - — r~  diag  (irsr*  >  uV~  »  ...  *  u~"~tt ~) 


"Ik 

\  1 


Vti  v  >  u  u  »  •  ♦  •  •  ii  u 
Nl"2  "lN3  VrF 


s3/n  h0n  t 

jg.i ,  u  UT 


-V  - B  X  -  U  IT 

lsa<S<«H  (Uj!?  aS'<  aBlr 


l” 


(31) 


#  a 


where  U  _  is  a  column  vector  whose  rows  are  labeled  with  the  same  scheme  as  is 
aBy 


used  for  the  elements  of  Dc>  and  whose  element  in  the  ij  row  is 


Then  the  matrix 


where 


U  „  (ij)  *  6  .  6,.  •  i  .  i  .  ♦  j 
aBy  v  ol  BO  oi  yO  Bi  yj 


L  J  FGF  du  =  K  -  l  H  U  l£  , 
B  l<a<B<y<M  aBy  uBy 


u  _  r  2  S7N  HflN 
aBy  '  B  u  ~ — a  du. 

NiP 


(32) 


(33) 


(34) 


4.9  Recursively  applying  the  matrix  inverse  lemma  [6]  to  the  right  side  of 
equation  (33),  the  inverse  can  be  written  as 


(/«  <02  F6F  do;)"1  =  K'1  +  l  K"1  U  -  M  n 
8  lso<3<y<H  aBy  aBy 


(35) 


where  for  the  purposes  of  this  paper  it  is  not  necessary  to  specify  further  the 

matrix  M  .  .  From  the  relations  defining  the  matrix  A  and  the  vector  U  „  ,  it 
aBy  aBy 

follows  that 


A'  U  .  =0. 
aBy  — 


(36) 
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4.10  Thus,  the  Inverse  of  the  covariance  matrix  for  the  Gauss-Markov 
estimate  for  D,  with  the  correlator  inputs  filtered  according  to  equation  (28),  is 

aT  pt  A  =  h  aT  k  (/b  “2  FGF  d“)_1  ** 

*  5T  *T  W  (37) 

-  (FIM). 

4.11  Thus,  the  correlator  system  optimally  filtered  according  to  equation 
(28)  provides  an  efficient  estimate. 
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Chapter  5 

Suboptimally  Filtered  Correlator  Systems 


t 


5.1  For  diverse  reasons  the  decision  not  to  use  the  optimal  filters  of 
equation  (28)  at  the  input  to  a  correlator  delay  measurement  may  be  made.  It  is 
then  relevant  to  investigate  the  degree  to  which  the  delay  estimate  is  degraded. 

The  question  can  be  answered  in  a  simple  way  under  the  following  hypothesis: 

a.  The  ratio  S/N-  is  the  same  at  each  sensor. 

b.  Identical  filters,  F,  are  used  on  each  input  channel. 

c.  The  suboptimally  filtered  Gauss-Markov  delay  estimate  is  used. 

Under  +^ese  hypotheses,  the  matrix  FSF  becomes 

F6F  =  | F | 4  (N2  +  HNS)  I 

W4  S"  %  ULr-  (38) 


l<a<B<YSH 


where  I  is  the  identity  matrix  and  UagY  is  defined  in  the  same  way  as  in  the 
preceding  chapter.  If  the  Gauss-Markov  estimate  is  formed  from  the  suboptimally 
filtered  delay  estimates,  the  inverse  of  the  covariance  matrix  for  the  Gauss-Markov 
estimate  is 

AT  PE 1  A  =  fc  K  (/b  “2  FGF  d“)_1  w 

=  (/B  "2  fp|4  (N2  +  "»)  d“)-1  aT  a  (39) 

=  T  f  S.|F|2,dM)2  aT  A  * 

**  JB  J2  |Fj4  (N2  +  MS)  du 
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The  FIM  for  this  case  is 

2  2 

(FIM)  =  ^  (/B  w2  dm)  AT  A.  (40) 

1  +  M^ 

5.2  Thus,  the  covariance  matrices  for  the  optimally  and  suboptimally 
filtered  estimates  differ  by  a  constant  factor,  and  it  is  easy  to  take  account  of 
the  effects  of  suboptimally  filtering  the  inputs.  Equations  (39)  and  (40)  can 
also  be  used  to  determine  the  degradation  of  the  delay  estimate  due  to  an 
imprecise  knowledge  of  either  S(ui)  or  N(u). 
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